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Abstract 

Selection systems and the corresponding replicator equations model the evolution of repli- 
cators with a high level of abstraction. In this paper we apply novel methods of analysis 
of selection systems to the replicator equations. To be suitable for the suggested algorithm 
the interaction matrix of the replicator equation should be transformed; in particular the 
standard singular value decomposition allows us to rewrite the replicator equation in a con- 
venient form. The original n-dimensional problem is reduced to the analysis of asymptotic 
behavior of the solutions to the so-called escort system, which in some important cases can 
be of significantly smaller dimension than the original system. The Newton diagram meth- 
ods are applied to study the asymptotic behavior of the solutions to the escort system, when 
interaction matrix has rank 1 or 2. A general replicator equation with the interaction ma- 
trix of rank 1 is fully analyzed; the conditions are provided when the asymptotic state is a 
polymorphic equilibrium. As an example of the system with the interaction matrix of rank 2 
we consider the problem from [Adams, M.R. and Sornborger, A.T., J Math Biol, 54:357-384, 
2007], for which we show, for arbitrary dimension of the system and under some suitable 
conditions, that generically one globally stable equilibrium exits on the 1-skeleton of the 
simplex. 

Keywords: replicator equation, selection system, singular value decomposition, Newton 
diagram, power asymptotes 

AMS (MOS) subject classification: 92B05, 92D15, 34C05, 34D05 

1 Introduction. Selection systems and replicator equations 

The evolution of replicators, which are the basic entities of the theory of natural selection, can 
be described with a high level of abstraction by means of the so-called selection systems or 



* Corresponding author: tel.: +1 (301) 451-6722; e-mail: karev@ncbi.nlm.nih.gov 
^e-mail: anovozhilov@gmail.com 
'e-mail: fberezovskaya@howard.edu 



1 



corresponding replicator equations, which, in their most generic form, can be written as follows. 

Let us suppose that every individual of a population is characterized by its own value of pa- 
rameter w, and uj takes its values from a measurable space (O, P) where is the set of admissible 
parameter values and P is a given measure. The parameter uj specifies an individual's invariant 
property and in most applications takes values in a finite set or in a domain of n-dimensional 
Euclidian space. Let us denote l{t,uj) the density of individuals with a given parameter value u 
with respect to the measure P; the distribution of this parameter can be continuous or discrete, 
depending on the nature of the problem we seek to describe with our mathematical model. 

An ab stract se l ection system (or, in the author's terms, a system with inheritance) was 



studied in iGorbanI (j2007l ) (see also references to earlier work therein) where a general selection 
theorem was proven. Roughly speaking, one of the statements of the theorem is that an infinite- 
dimensional abstract system with inheritance "tends" in the course of time to a finite dimensional 



system (see Gorban ( 200?! ) for the exact statements). This result justifies the special attention 



to selection systems with a discrete distribution of the parameter. The simplest example of 
a discrete distribution of parameter u is given by interpretation of u as merely an index of 
interacting subpopulations; in this case l(t,u}) is naturally to interpret as the size of the w-th 
subpopulation. To emphasize the discrete nature of the distribution of uj in some problems, we 
will use it as an index: luj{t) (or, more traditionally, li{t), replacing w with the index i). 

If we denote the per capita growth rate of w-th subpopulations as F{t,uj) and assume the 
overlapping generatioi is and smooth ness of l(t,uj) in t for each fixed u, we obtain the abstract 



selection system (e.g.. iGorbanI (|2007l )) 



^l{t,co) = l(t,Lo)F{t,co), l{t,co)>0, (1.1) 

where the initial condition /(0,a;) is given, and the growth rate, or fitness, F{t,u}) can depend, 
among other things, on the total population size N{t) = J^l{t,uj) du (where the integral is 
replaced with the sum if the distribution of uj is discrete). The exact form of F{t,uj) we will be 
working with is given below. 

It is straightforward to infer that the frequencies of subpopulations, 



N{t) ' 

satisfy the replicator equation: 

^^p{t,uj)=p{t,io){F{t,uj)-Et[F]), (1.2) 

where = j^F{t,uj)p{t,uj) du denotes the mean fitness of the total population at the 

time t. The natural phase space of the replicator equation ()1.2p is given by {pit^u): p{t,uj) > 
0, f^p{t,uj) duj = 1}, in the discrete case it is the simplex 5„ = {piif) : Pi{t) > 0, ^^iPiit) = 1} 
(here and below we assume that generally there are n interacting subpopulations, and the 
notation means Yl'i=i)- remark that the same replicator equation (|1.2p can be obtained 
for different selection systems (jl.ip . e.g., it is true if the growth rates in two selection systems 
(jl.ip differ by a function that depends only on the total population size; when passing in the 
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opposite direction, from the replicator equation to the selection system we always choose the 
simplest one. 

Naturally, equation ()1.2p should be supplemented with the equation 



dt 



Nit) = Et[F]Nit), 



if the fitness F(t,u}) depends on N{t). 

The replicator equation (|1.2p comprises well-establish e d biq r nathe n iatical models in quite 
distinct evolutionary contexts, see lHofbauer and Sigmundl (|l998l . hmi ): ISchuster and Sigmundl 
(|l983l l. We survey these models briefly. 

One of the first replicator equations was used, at least implicitly, by Ronald Fisher, John 
Haldane, and Sewall Wright to study the evolution of multiallelic one-locus gene frequencies 
under the forc e of natural selection in a sexually reproducing diploid population (for more 
information see lHofbauer and Sigmund (|l998l . l2003l ll. If a gene locus with n alleles is considered 
the frequency of the i-th allele is denoted as pi, and the Hardy- Weinberg equilibrium is assumed, 
then, in the usual way, the replicator equation is obtained: 



5«« 



3i 1^ 



(1.3) 



where rriij is the Malthusian fitness of genotype with alleles i and j. Given the assumptions, we 
have that in this case m. 



rrii. 



, the fitness matrix M = {rriij} is symmetric. In equation ()1.3p 
the intrinsic growth rate of the i-th allele depends linear on the frequencies of other alleles, in 
our notations, 

Fuj{t) = ^ . m^jpj{t), LO = 1, 



If the fitnesses rr iij ar e const ant, then it is known that 



(|l998l ): IShahshahanl (|l979l 'l: ISvirezhev and Pasekovl ( 
lation is monotonically increasing. 

Sometimes it is natural to assume that the corresponding growth rates depend also on the 
population size and therefore include the density-depe ndent effects in the model, e.g. , similar 



is a gradient system ([Hofbauer and Sigmund 
1990)), and the mean fitness of the popu- 



to the well-known Verhulst-Pearl logistic equation (e.g.. ICharlesworthI (ll97lh:lGinzburg| (1197711). 
One p articular model with explicit birth and death terms, considered in lDesharnais and Costantino 
(|l983l ). has the form 



hi,f{N)-di,g{N), 



:i.4) 



where bij and dij are the per capita density-independent rates of recruitment and mortality, 
respectively, associated with the genotype {i, j}. Therefore to describe the evolution of the 
allele frequencies it is necessary to consider the following problem 



d_ 

lit 
d_ 

lit 



= Pi {{hi - Et[h])f{N) + {Et[d] - di)g{N)) 
N = N {Et[h]f{N) - Et[d]g{N)) , 



1, 



:i.5) 



where bi = Y^jhjPj, di = YjjdijPji Et\h] = Y^^hpi, Et[d] = Y^idiPi, and the explicit depen- 
dence on t for Pi and N was suppressed for simplicity. System ()1.5p can be replaced with the 
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following selection system 



-h = k (bifiN) - dig{N)) 
at 



,n, 



N{t) = Y,.h{t)- 



Another example of the replicator equation is given by the equation for the frequencies of 
the pure strategies in the population w it hin the framework of the evolutionary game theory 
fsee iHofbauer and Sigmi^ (|l998l . hooA ): IXavlor and .Tonke^ ^197^ )). If it is assumed that Uij 
denotes the payoff of the player with pure strategy i against the player with pure strategy j, 
then the dynamics of the frequencies of the players in the population are given by 



^Pi(t) = Pi{t) aijPj{t) - ^ ajkPj(t)pk{t) 



,n, 



:i.6) 



where, contrast to (|1.3p . the matrix A is an arbitrary re a.1 n xn matrix . In the case of a contin- 
uum o f pure strategies an analogue to ()1.2p is obtained (Bomze ( 199d ): Hofbauer and Sigmund 
( 2Q03I )). System ()1.3p is a particular case of ()1.6p with a symmetric A, such matrices describe 
partnership games. The selection system, corresponding to (jl.6p . has the form 



hit) (Y^^ 



Another well studied particular case of ()1.6p is the so-called hypercycle equation, which 
is given by setting Ojj = ki \i j = i — 1 and aij = otherwise. The hypercycle equation de- 
scribe s a catalytic loop of self-re plicating macromolecules, each promoting replication of another 
type (jEigen and Shuster Completing the list of application of the replicator equation 

(|l.'2p we note that the classical equations of Volterra describing dynamics of interacting pop- 
ulations can be transformed into the form ()1.2p by means of an invertible change of variables 



(see iHofbauer and Sigmundl (|l998l )). However, we remark that the methods, described in the 
present text, are better to apply directly to the Volterra systems. 



Ef fective methods of analysis of selection systems (jl.ip were developed recently (see iKarev 
1,2010') and references therein) for particular form of the fitness function F{t, cj). Here we present 
the applications of these methods to the selection system (jl.ip and, consequently, to the repli- 
cator equation (|1.2p . It turns out that some of the systems obeying the replicator equations can 
be effectively analyzed and solved even for large n. 

Our paper organized as follows. In the next section we present an algorithm, which allows us 
to replace a given selection system with an equivalent problem. This equivalent problem in some 
particular cases can be of significantly lower dimension than the original system, and this is the 
case when the suggested methods should be taken advantage of. Section 3 is devoted to the 
methods how to transform a given replicator equation so that the methods of Section 2 can be 
directly applicable. In Section 4 we present a non-trivial application of the proposed technique 
to the replicator equation, which is supposed to model the evolution of sensory systems, and 
obtain a general proof, under some suitabl e conditions, of a conjec t ure, w hich was proved only 
for particular cases in the original study lAdams and Sornborgerl ()2007l ). The last section is 
devoted to conclusions, and Appendix contains some auxiliary facts. 
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2 How to solve selection systems 



Here we present an algorithm that ahows us to formaUy write down an expHcit solution to the 
selection system (jl.ip . which we rewrite here for convenience, 



d 

— l{t,uj) = l{t,uj)F{t,uj), l{t,uj)>0, 



when the fitness F{t, u) has the following special form: 

'r- 

where Gj, Kj are the so-called regulators 



Ui{t, Gi)ipi{io) + 2^ .^^ Vj{t, Hj)Tpj{uj), 



(2.1) 



Gi{t) = / gi{uj)l{t,uj)dw = N{t)Et[gi], i = l,...,mi, 

Hj{t)= / hj{io)p{t,u))dio = Et[hi], j = 1, . . . ,m2, 
Jn 



(2.2) 



Ui, Vj, Qi, hj, ipi and ipj are given functions, mi, m2 > are constants, andp{t, to) = l{t, u))/N(t). 
The probability density function p{t,uj) solves the replicator equation ()1.2p . 

In applications functions (pi{u}), can be interpreted, for instance, as particular pheno- 

type traits that characterize an individual with the parameter value w; Ui{t,Gi) and Vj{t,Hj) 
then describe the contribution of the corresponding phenotype traits to the fitness (mean num- 
ber of descendants per individual) at the time moment t provided the values of regulators Gi and 
Hj. Note that the traits {(pi} correspond to the density-dependent regulators Gi, whereas the 
traits {V^j} correspond to the frequency-dependent regulators Hj. We divide the regulators into 
these two group for convenience, although it should be clear that the theory could be written 
only for Gi. 

This particular form of the s election system ()l.ip . ()2.ip com prises niany meani n gful mathe- 
matic al models (s ee^Karev^ ('20ld) for the general theory and, e.g., Karev ( 2003 . 20051 ): Karev et al 
(|200fil ): lNovozhilov (.2004 . 200S ) for various applications). Model ([11]) , ([21]) , ([22]) defines, in gen- 
eral, a complex transformation of the initial distribution l{0,u}). The remarkable fact is that 
model (|l.ip . (|2.ip can be reduced to an equivalent system of ordinary differential equations 
(ODEs). Here we present only the algorithm; the proofs can be found in iKarevI (|2O10l ^. 

We introduce the functional on the space of measurable functions of the parameter co: 



M{z; \,5) = j^ z{u:) exp {Y^^^ + Y^li '^iV'i(w)} p(0, 

= Eo ^zexp |2^.^^ \m + l^-^^ '^iV'j|J , 

where p(0,w) = /(O, a;)/iV(0), and A = (Ai, . . . , AmJ, ^ = (Si, . . . ,6m2)- 

Next we write down the escort system of ordinary differential equations: 



uj) duj 



(2.3) 



d_ 

di 

d 

di 



i(t) =u,(t,iV(0)M(5,;q(t), s(t))), gi(0) = 0, i = l,...,mi, 
Sj{t)=Vj{t,M{hj;q{t), s{t))/M{l; q(t), s(t))), Sj{0) = 0, j = l,...,m2, 



(2.4) 
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where q(t) = . . • ,gmi(t)) and s{t) = {si{t), . . . ,5^2(0)- 

Using the functional (|2.3p and the solutions to (j2.4p we can write the solution to the selection 
system ([Ll]) , ([21]) as 

l{t,co) = l{0,Lo)K{t,Lo), (2.5) 

where 

E:(t,a;) =exp|J_;^^^g,(t)(^,(a;) + 2_;^.^^s,(t)Vi(a;)|. (2.6) 
We also have that the total population size is equal to 

N{t) = iV(0)M(l; q(t), s(t)); (2.7) 
the values of regulators are given by 



G,{t) = N{0)M{gi; q(t), s{t)), i = l,...,m,, 

_ Mjhy, cijt), sit)) {2.i 
- M(l; q(t), s(t)) ' ^-l'---'-^' 



and the current probability density function of the parameter distribution can be presented in 
the explicit form as 

Formula ()2.9p is the central result of the theory; it gives the solution to the replicator 
equation p.2p and allows us to compute all the statistical characteristics of the underlying 
parameter distribution in the self-regulated selection systems of the form (jl.ip , ()2.ip . 

The major technical tool in the considered approach is the functional M{z; A, S), which is 
formally well-defined for any given initial distribution p(0,u}); in practice, however, it might 
be difficult to evaluate M(z; A, S) for particular functions z{uj). A significant simplification is 
achieved if fitness F{t,uj) depends only on the regulators of the following form: N(t), Et[fi], 
or N{t)Et[(pi]. In this case it is straightforward to see that instead of the general functional 
M{z; A, S) we can use the moment generating function (mgf) of the initial distribution p(0, u), 
which is defined as 

Mo{5) = Eo exp|^^'^(^i(/7i} . 

Indeed, 

M(l; q(t)) = Eo [exp {Y^^i = ^'^^Ht)), 

M((/?fc;q(t)) = Eo [(^fcexp{^^'^(Zi(t)^i}] = _Mo(q(t)). 

The same holds for the frequency dependent regulators. 

Using (|2.10p the right hand side of ()2.4p can be rewritten in terms of the moment generating 
function of the initial distribution, which is supposed to be given. Remark that the moment 
generating functions are known for many important probability density functions. 

As a simple example we consider the following 
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Example 1. Let the Malthusian fitnesses be multiplicative, i.e., ruij = mirrij for the given 
vector m = (mi, . . . , m„), where rrii 7^ nrij for any i 7^ j. The asymptotic outcome of the 
dynamics is well known and trivial: all alleles but the one with the highest fitness rrimax Sire lost; 
here, additionally to this known result, we obtain also a simple expression that can be used to 
compute time-dependent behavior. System ()1.3p for the multiplicative fitness can be rewritten 
as 

-Pi{t) = Pi {m,Et[m] - {Et[in]f) , i = l,...,n, (2.11) 

where -Bt[m] = YlJ=i nT'jPji't)- The following selection system corresponds to replicator equation 
(|2TT]) : 

^k{t) = k{t)miEt[in\, i = l,...,n, (2.12) 

and belongs to the class (jl.ip . ()2.ip . 

As before, denote Mo(A) the moment generating function of the initial distribution Pi{0), 
Mo(A) = exp{Amj}pj(0) = £'o[exp{Am}]. The escort system ()2.4p consists only of one 
equation: 

The solution for the frequencies is given, according to ()2.9p . by 

p.(t)=p,(0) "!P|"!;;^y , Eo[K{t,-)]=Y.-^Ms{t)m,}pl (2.14) 

where p^ denote the initial conditions, p^ = Pi{0) > for any i. 

Prom (j2.13p . and using the change of the variable s{t) = — Inn(t), we obtain the equation 
for the new variable u 

■ _ ^i'lTlip'^U^* 

where di = rrimax — rrii, which implies that one of di = 0. From the last equation it follows that 
u{t) — >■ as i — >■ 00, which yields that s{t) — > 00. Using the last fact and the solution to ()2.13p 
we have 

— — = -fy exp{s(t)(mi - mj)\ 00 
Pj{t) p) 

if rrii > rrij, which is possible only ii pj{t) due to the constraint YliPii^) — ^■ 

The major advantage of the considered approach is that if one needs the time-dependent 
behavior of system ()2.1ip then, instead of solving n differential equation it is suffice to solve 
only one differential equation (|2.13p for the auxiliary variable s{t). In a similar vein the case of 
additive fitness rriij = mi + nij can be analyzed. 

3 Reduction of a general replicator equation by means of matrix 
decompositions 

It is obvious that only in exceptional for the system ()2.1ip , the theory of Section [2] can 

be applied to the replicator equation ()1.6p directly. In the general case we need to rewrite the 
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interaction matrix A from the equation (|1.6p in the form, suitable for the described technique. 
In this section we propose a method to apply general technique of the analysis of the selection 
system (jl.ip to the replicator equation. 
We start with a symmetric matrix A. 



3.1 Spectral decomposition 

Let us consider again the equation for the allele frequencies in diploid population ()1.3p with a 
constant matrix M. According to the interpretation of this equation, matrix M is symmet ric, 
niij = niji. Any symmetric real matrix M can be presented in the form (e.g., lOrtegal (jl987l )) 

M = Aihih[ + A2h2h5 + . . . + XkKhl, (3.1) 

where A,, i = 1, . . . ,k are the real eigenvalues of M, k is the rank of M, hj, i = 1, . . . ,k are 
the corresponding right eigenvectors that satisfy h^hj = 1, h^hj = 0, i ^ j, and r denotes 
transposition. The form ()3.ip is the spectral decomposition of M. If we denote the j-th element 
of the i-th eigenvector hj as hji then each element of M has the form 

niij = Xihiihji + \2hi2hj2 + . . . + Xkhikhjk- 

According to the last equality, system ()1.3p takes the form 

j^Piit) = Piit) (Y'^^^i >^Jh^JEt[h,] - Y".^, XjiEtii^j]?^ , i = 1, . . . , n, (3.2) 

and hence belongs to the class of selection systems (jl.ip with fitness ()2.ip . 
Consider the mgf of the initial distribution of h^: 

The last expression yields the following escort system: 

±^ ^ ^ Eo[hj explX]^^^ Sm{t)h^}] 
dt'^' ' Eo[eME'm=iSm{t)h^}] 

_ ^ . Yli=lPi^ij J2m=l ■^m{t)him j = \ ^ 



System p.3p can be rewritten in the compact form 

d , d_ 

Hence the solution to system ()1.3p is 



(3.3) 



^Si = Ai^lnMo(s), j = l,...,k. (3.4) 



P.(i)=P.(0);g^^, K(t,i)=exp|j]'^^s,(t)/i,,|, Eo[K{t,-)] = MoHt)). 
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To sum it up, we replace n-dimensional problem ()1.3p with A:-dimensional system (|3.3p . We 
point out that the latter system can be no easier to solve than the original one. The suggested 
approach is beneficial only when A: <C n (here we only speak of finite dimensional systems (jl.Sp ). 

The spectral decomposition approach as in the previous example allows us to reduce the 
original problem ()1.6p to the system in the form (jl.ip . ()2.ip only for symmetric A. In the 
gener al case with an arbitrary real matrix A we can apply the singular value decomposition 



(e.g., iJoMd (12002)). 

3.2 Singular value decomposition 

Well known that, given an arbitrary matrix A of dimension n x n, A can be written as 

A = USX^, (3.5) 

where U, X are nx k matrices, each of which has orthonormal columns so that U^U = X'^X = 
Ifc, where Ifc is the identity k x k matrix; H is a k x k diagonal matrix with non-n egative 



eleme nts cxi > cj2 > . . . > cr^ > on the main diagonal; and k is the rank of A (see, e.g., iJolliffe 



(|2003)). The representation ()3.5p is called singular value decomposition (SVD). Singular values 



aj, i = 1, . . . ,k are the square roots of the eigenvalues of AA^ (or A^A), which are given by 
(rf > a2 > ■ ■ ■ > > 0; columns of U and X are the right eigenvectors of AA'^ and A'^A 
respectively; each column corresponds to its own (t|. 



Prom (|3.5p it follows that 

m=l ■' 

where Uim^ are the elements of U and X"^ respectively. Using this representation, any matrix 
A in ()1.6p can be written such that the fitness F{t,uj) in (|1.2p is in particular form (|2.ip . That 
is, we obtain 

^K(i) = Pi{t) i X]j=i UijCrjEtixj] - ^.^-^ ^^jEtiujjEtlyij] j , (3.6) 

where Uj, Xj denote the j-th columns of matrices U and X"^ respectively. 
Singular value decomposition enables us to split elements aij into parts 

'^ira^mXjrm ^ — 1, . . . , fe. (3.7) 

If only q < k such parts are retained, then the expression 

El 
171=1 ■' 

provides an approximation to aij in a sens e that na,;,- gives the best possible rank q approximation 
to aij (the proof can be found in iGabriell (|l978l l) when 



"J 

is minimized with respect to any matrix (qaij) of rank k. 
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The last point bears a close relationship with the principal component analysis (jjolliffe 



( 20021 )). Let us assume that matrix A represents a table of n observations of m random variables. 



The usual way to reduce the dimensionality of this data-set, while retaining as much as possible 
variation presented in it, is to apply the principal component analysis. It can be shown that 
parts ()3.7p give the contribution of the corresponding principal components to the data-set; 
retaining only q parts corresponds to retaining q first principal components. 

We remark that the singular value decomposition essent ially depends on the choice of the 
scalar product (see, e.g., Ch. 5 in Gorban and Karlin ( 20051 )). In the case of high- or infinite- 



dimension models we can, taking a finite number of components, reduce the original model to a 
finite dimensional selection system with the fitness in the particular form (j2.ip . In what follows 
we study only finite dimensional models of the form (6) with the standard scalar product. Note 
that if matrix A has the rank k, then the escort system (j2.4p is A:-dimensional. 

It is reasonable to assume that matrix A in applications is known only approximately. Re- 
taining q first principle components would correspond to the approximation of the matrix A 
with the best possible matrix of rank q; from the standpoint of the theory presented in Section 
[21 such approximation is useful because the dimension of the escort system is reduced to q. 

Let us rewrite system ()1.6p in the form 

^K=p,((Ap),-p^Ap). (3.8) 

The stationary points are found as the solutions of 

Ap = (p^Ap)l„, 

where 1„ = (1,1,...,1)'^. Using ()3.5p the last system can be rewritten as 

ai^t[xi]ui + . . . + (Tfc^t[xfc]ufc - (p^Ap)l„ = 0, (3.9) 

which shows that in general an isolated polymorphic equilibrium can exists only when k = n, 
or when vectors ui , . . . , Ufc, I n are linear dependent (cf. with the proof of similar conjecture in 



Sornborger and Adamsl ()2008l )). On the other hand, when /c < n it is possible to have manifolds 
of non-isolated equilibria for which -Et[xi] = . . . = -Et[xfe] = (see Section [33] for a particular 
example). Therefore, approximation of the original system with a matrix of smaller rank can 
yield either lost of information on isolated polymorphic equilibria, or appearance of manifolds 
of non-isolated equilibria. 

On the other hand, it is reasonable to expect that in applications matrix A can have exactly 
rank k, whereas its estimate A, which is known to the researcher, can have the rank n due 
to, e.g., noise effects. In this case the reduction technique does not loose any information, 
and, additionally, discard the fictional information, which appears in the matrix A thanks to 
the estimate errors. Again, using analogies with the principal component analysis, we can 
use known technique to infer the dimension (the rank of the matrix) that contains the principal 
informat ion (for a review article how to determine the number of significant principal components 
see, e.g.. ICangelosi and Gorielv (|2003)). 



It is interesting to note that for the case of a symmetric A the escort system is a gradient 
system. 
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Example 2 (Partnership games and gradients). In Example 13.11 we already considered the 
replicator equation with symmetric matrix using the spectral decomposition. Here we apply 
SVD to such systems. 

Let = A, then (|3.5|) becomes A = USU'^ since AA'^ = A'^A, and each element of A 
has the form 

m=l ■' 

Denote where Uj is the i-th column of U. Then we have that A = UU , and 

aij = X]m=i ^imUjm- Finally, the replicator equation takes the form 



L, . . . ,n, 



with the corresponding selection system 

^k{t) = k{t)Y,''.^^UijEt[uj], i = l,...,n, (3.10) 

which is as required by (|2.ip . As before denote Mq{S) the moment generating function of the 
initial distribution of the elements of uj. The escort system now reads 



d d_ 
or simply 



= ^lnMo(s), Sj{0) = 0, j = l,...,n, (3.11) 



s = VlnMo(s), s(0) = 0, 

which is a gradient system with the potential — InMo(s) in the usual space with the standard 
Euclidian metric. 

3.3 Analysis of the replicator equation having the matrix of rank 1 

To conclude this section we consider the problem of finding all possible asymptotic states of 
the replicator equation with the interaction matrix that has rank 1 (note that Example [T] is a 
particular case of this problem). 

According to SVD any matrix of rank 1 can be presented as A = ab"^, where a and b are 
vectors. 

It is a simple matter to determine the asymptotic states in the replicator system with such 
matrix A in the case b > 0. 
The escort system reads 

Y^ibiPi e^p{ais} 



EiP°exp{ais} 

where all > 0, and the initial condition s(0) = 0. 

Using the change of the variable u = exp{— s} we obtain that 



(3.12) 



(3.13) 
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where di = maxj{aj} — Oj. Note that at least one di = 0, so that the previous equation can be 
rewritten as 

u = -u^^ 0^ d ' 3.14 

where Ik = {i- ai = ma,Xj{aj}}, and all dj > if i ^ /fc- From ()3.14p and if b > it follows that 
the origin is a grobally asymptotically stable equilibrium and, when t ^ oo 

«(t)^exp{-At}, A = %^^, 

Z^iG/fe Pi 



which means that s{t)/t \ as t ^ oo. 
We have that 



Piit) 



p° expjaj s (t) } p'ju'^^ 



which yields that 



Pj{t) ^ ^ ' if^e/fc, (3.15) 

and Pj{t) —7- otherwise, which is generalization of the result in Example [TJ Note that if b < 
a similar result is valid if we denote 1^ = {i- a-i = minjjaj}}. 

To analyze the general case, when elements of b can have arbitrary signs, we introduce two 
parameters: = Y^ieh ^^Pi' ^2 = Ei = Eo[h]. 

If ^1 > and ^2 > then, due to continuity of the right hand side of ()3.14p . equation ()3.14p 
may possess zero or even number of equilibria {ufc}- Using the fact that n = is a stable 
equilibrium and recalling that the initial condition is n(0) = 1, we obtain that u{t) — ?> or 
u{t) — 7> u* when t oo, where u* is the closest equilibrium of (jS.lSp to u = 1 belonging to 
(0,1). 

In the case < Oi ^2 > there is always equilibrium u* € (0,1), and -u = is unstable, 
which means that u{t) u* when t oo. 

In the case > 0, .^2 < it is possible that u{t) 00 when t ^ 00 if there are no equilibria 
of (|3.14p when u € (1, 00), or u{t) — >■ U* , where U* is the closest equilibrium of (|3.14p belonging 
to (1, 00). The case ^1 < 0, ^2 < is analogous to the previous one. 

Therefore, we showed that u{t) can tend to 0, w or cxd when t 00. Using this fact, the 
explicit expression for frequencies, and auxiliary notation J = {j : aj = minjjaj}}, we obtain 

Proposition 1. Three types of asymptotic behavior of the solutions pi{t) to the replicator equa- 
tion (|1.6p are possible, if the interaction matrix of the replicator equation has rank 1: 
1) If u{t) — )• as t —)• 00 then 

Pj(^) ^ — ^ 

l^i&Ik Pi 

and pj(t) — )• otherwise; 
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2) If u{t) — )• n as t ^ oo then 



3) If u{t) — > oo as t oo then 

l^iejPi 

and pj(t) — )• otherwise. 

In Proposition [1] we studied the generic case when Eq [b] ^ 0. Note that all distributions p 
such that Eq [h] = are equilibria of the replicator equation 

j^Piit) = p^{t){aiEt[h] -Et[a]Et[h]) i = l,...,n. 

Therefore, the (n — 2)-dimensional subset of the simplex Sn, = {p: p G Sn, EQ[h] = 0} 
consists of interior equilibria of the replicator equation. It is worth pointing out that if u{t) — >■ 
u 7^ then the limit distribution p(oo) is an equilibrium belonging to S^. Indeed, in this 
case s{t) — > exp{— -u} < oo, while the variable s{t) was defined by the equation s = Et[h], 
s{t) = E-rlh] dr. Hence, s{t) is bounded for all t only if Et[h] — )■ as t — )• cxd and -EoM = 
for the limit distribution p. 

Proposition [1] not only shows that the asymptotic states of the replicator equation with the 
matrix of rank 1 can be only equilibria, it also points out that the case of polymorphic (interior) 
attracting equilibrium, albeit non-isolated, is not exceptional for the general vectors a and b. In 
particular, if the parameters defined above are such that < and ^2 > then the asymptotic 
state is always polymorphic. 



4 Analysis of a class of replicator equations 

In this section we consider a non-trivial example, where the application of the suggested methods 
allows us to give a proof for a problem concerning the evolution of sensory systems. 

Example 3. Motivated by a problem in the evolution of se nsory systems where gai r is obt ained 



by improvements in detection are offset by increased costs, lAdams and Sornborgeii (|2007l ) con- 



sidered the dynamics of the replicator equations ()1.6p with the matrix of the form 

aij = aibj - Ci, i,j = l,...,n, (4.1) 

where a = (oi, . . . , UnY, b = (61, . . . , bnY, and c = (ci, . . . , CnY are given non-negative vectors. 
They showed, using topological arguments, that in the case of general position (see below) and 
for dimension n < 5, the system can have only one global attractor, and this global attractor is an 
equilibrium having at most two non-zero components. We prove this conjecture, using completely 
different methods, for an arbitrary n, with some additional (mainly technical) conditions on 
a, b, c. 
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Using the notations from Sections [T] and [2l we obtain the rephcator equation in the form 



j^Mt)=piit)iaiEt[h]-Ci-Et[F]), i = l,...,n, (4.2) 

and the corresponding selection system is 

^k{t) = k{t){aiEt[h]-Ci), i = l,...,n, (4.3) 

to which the methods from Section [2] can be applied. To clarify the connection with (|2.1|) we 
write down explicitly 

V'i(w) = a(a;), vi{t,Hi) = Hi, Hi = Et[h], 

where to takes the values from the discrete set = {1, . . . ,n}. Therefore, the escort system 
()2.4p has the form 

^ „ _ Eo[hexp{asi{t) + cs2{t)}] n 
dt''^'>~ Eo[eM^si{t) + cs2{t)}] ' "1^^^-^' 

j^S2{t) = -1, S2(0) = 0, 

or, integrating the second equation in the last system and dropping index of si{t), finally we 
obtain one differential equation 

d _ ^o[bexp{ag(t) - ct}] _ Y,7=i exp{ajg(t) - at} (r^^^r. 

dt ^' EoieMMt) - ct}] Er=iP°exp{a,s(t)-c,t} ' ^ ^ ' ^ ' ^ 

where all > 0. 

Before analyzing equation ()4.4p we note that we consider only the generic case for the game 
matrix given by (|4.1|) . The genericity condition in our case reads as follows: any projections 
of the three vectors a, c, and In to any of the three dimensional subspaces of spanned by 
three standard coordinate vectors are linearly independent (this means that isolated equ i libria of 
the system can have at most two non-zero coordinates, see also lAdams and Sornborger jioo^)). 



Putting in other words, this condition means than for any indexes i,j, k the following holds: 




or 

«i(cj - Cfc) + aj{ck - Ci) + afc(cj - Cj) / 0. (4.5) 

We are particularly interested in the limiting behavior of Pi{t) as t — > oo. Here we show, 
analyzing the escort system of our replicator equation, that Pi{t) tend to an equilibrium of the 
initial replicator equation as t — )• oo and that this equilibrium is the global attractor of our 
dynamical system for arbitrary n. 
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First, we make the change of the variables 



] I (4.6) 
V = expj— t} ^ t = — Inu. 



In the new variables equation ()4.4p takes the form 



(4.7) 



where di = maxjjaj} — a^. Note that at least for one i di = 0. We can also assume, without loss 
of generality, that minj{cj} = (in general, we scale q = Cj — minjjcj} and drop the tilde for 
notational simplicity) . 

The initial condition for (j4.7p is u{v = 1) = 1. 

We rewrite (|4.7p as a dynamical system on the plane 



u = u 



^ — '1=1 



n(0) = t;(0) = 1, 

where the derivatives are taken with respect to some dummy "time" variable. We remark that 
system ()4.8p has an isolated equilibrium 0(0,0), and the axes u = and v = are orbits so 
that 0(0,0) cannot be monodromic (focus or center). Using the function / = (u^ + f^)/2 we 
find that 

Ltf = ^u+^ij = u^y hip%'^^v^^ + y p%'^^v^^ > 

du OV 

for li > 0, w > 0. Here -Lf(-) is the derivative along the orbits of the dynamical system ()4.8p . 
The last expression implies that the first quadrant of the phase plane (the one we are actually 
interested in) o f (I4.8p is a repelling parabolic sector (i.e., it is a parabolic sector, see, e.g., 



Andronov et al.l (jl973l ) for the definitions, for which the origin attracts the orbits when "time" 



tends to — oo). 

The qualitative the ory of ordinary differen tial equations on the plane is an extensively re- 
searched area (see, e.g., Dumortier et al. ( 20061 )). especially in the case when the right hand sides 



are analytic functions. In our case this would mean that di and q are natural numbers, but it 
is a straightforward procedure for our problem to extend the basic necessary results to the case 
when Cj , S Q+ (see Appendix) , so in the following it is assumed that q , di are non- negative 
rational numbers. 

If point 0(0,0) is not monodromic therefore there are c haracteristic dire c tions along which 



the orbits of the dynamical system approach the equilibrium I Andronov et al.l (jl973l ). According 
to ()4.6p . the behavior of the orbits of ()4.8p when u,v ^ 0(0,0) determines the asymptotic 
behavior of s{t) when t — )• oo. In the following we shall call the orbits of ()4.6p 0-orbits if 
u,v ^ 0(0,0) for positive or negative "time" directions. We also recall that 0-orbits of system 
(j4.8p have a power asymptote with a positive exponent p and a non-zero coefficient C if 

n = Ot;^(l + o(l)), O / 0, /> > 0, V ^ 0. 
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Figure 1: The Newton diagram is given by the convex polygonal line passing through the points 
Aij^, . . . , Ai^ (an illustrative example). The diagram consists of 4 vertexes with the four indexes 
of the second kind /3j and of three line segments with the indexes of the first kind aj. See text 
for details 



The theory of th e power asymptotes o f 0-orbits of an isolated equ i librium in t h e plane is well- 

developed (see, e.g., Berezovskava ( 19761 ): Berezovskava and Kreitsei ( 1975 . 197G ): Berezovskava et al, 
(l2007l ): lBriu^(|l989l ll. and we apply it here following mainlv lBerezovskava and Kreitseil (j 19761 ). 

To summarize the main results we need the notion of the Newton polygon. We assume that 
all > 0. Introduce rectangular coordinates in the plane and to every bi ^ assign a point Ai 
with coordinates (ci,di). Consider the convex polygonal line J\f passing through points of the 
set {Ai} joining {0,dij) and (cij.,0) such that each Ai lies above or on J\f. This line is known as 
Newton's polygon or Newton's diagram (see Fig. [1]). The Newton polygon consists of a finite 
number of line segments Mj, whose angles with x-axis are between and vr/2, with end vertexes 
Ai- and Aj.^j, j = 1, . . . , K — 1, and K is the number of vertexes. To each vertex Ai. of the 
Newton polygon M we assign the index of the second type /3j = bi- , and to each line segment 
J\fj we assign the index of the first type 



(di, -di,+i) {ai^^,-ai^) 



Using the results from Berezovskava and Kreitser ( 19761 ) we have the following theorem: 



and di ^ dj, Ci ^ Cj for any i ^ j, bi > 0, > 



Theorem 1. Let di, Ci € Q+ , 6^ € II 

for any i in system (|4.8p . and bi^ 7^ ^ij+i f^''" "^2/ ^''^^ vertexes of line segments of the Newton 
diagram. Let condition (|4.5p hold, and index of any vertex Ai^ be not equal to the indexes of 
the adjacent edges, i.e., f3j ^ aj and f3j ^ a^+i. Then all O-orbits of system (|4.8p have power 
asymptotes. The positive exponents p of the power asymptotes of O-orbits can be 

i) p = /3j = bi- , where f3j is the index of the second type of the vertexes of the Newton polygon, 
if for the vertex Ai. aj < (3j < a^+i holds; 
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ii) p = /3i i//3i < ai; 

iii) p = /3k if Pk> aK^i; 

iv) p = aj, where aj is the index of the first type of the line segment Afj of the Newton 
polygon, if the function ^{z) = Pi.ihj —Oij)z '-^ +P^i^^^{bij+i —aj)z 'j+i has non-zero root z = C, 
and this root is the coefficient in the power asymptote. 



Remark. The major assumptions in Berezovskava and Kreitsei ( 19761 ) are different from those 



given in Theorem [T] but follow from them. In particular, for any line segment Mi it is necessary 
to consider two functions: 

n / \ di . Ci . +1 , di.,, Ci . , . +1 

Pj{u,V) = PijU 'IV J + ^ + ^V J + 1 , 

Q,{u, v) = bi^p^u'^'^^V'^ + hi^^,plu''^^+^+\''^^+^ , j = l,...,K. 



Then the conditions from iBerezovskava and Kreitser (jl97fil ^. adapted to our problem (|4.8p . 



can 



be stated as: the functions Pj{u, 1) and Qj(u, 1) cannot have common non-zero real roots, which 
follows from the fact that bi. ^ 6,^.^^ for any two vertexes of the line segments of the Newton 
diagram; and the function ^j{u) = —ajuPj{u, 1) + Qj{u, 1) cannot have multiple non-zero real 
roots, which follows from the fact that ^j{u) is a binomial, and should be identically zero to 
have multiple roots. The condition for Pj, and Qj to be binomials follows from ()4.5p . 

Noting that we are given the initial conditions u{0) = v{0) = 1 we conclude that there is only 
one orbit passing through the point (1, 1). This orbit has a power asymptote u = Cv^^l + o(l)) 
when V ^ (first quadrant is a parabolic sector, where all the orbits have power asymptotes), 
and the exponent p can be either the index of the first type or of the second type of the 
corresponding Newton polygon Af. 

Having the power asymptote we can rewrite ()3.2p in the variables u, v: 



„ - py'"""' (49) 



or, using u = Cv^{l + o(l)). 



For the following we need ( Berezovskava and Kreitser! ( 19761 )) 



Lemma 1. For any index of the second type j3j of the Newton polygon M , which can he an 
exponent in the power asymptote (see TheoremU\), we have 



di^ + PjCij = e| > 0, 
dk + f3jCk = e^j +el, k^i^, 



2,-2 , . . (4-11) 



where if, > 0. 
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For any index of the first type aj of the Newton polygon J\f we have 

dij + ajCi. = > 0, 
+ OjCi^+i = e] > 0, (4.12) 
dk + ajCk = e] + el, k^ij,k^ij+i, 

where el > 0. 

First suppose, without loss of generality, that p = (3j = hi. Then, from (j4.10p and Lemma 
[H we have 

Mv) = rf^!Hl±^(l» . ,4.13) 

p;c* (1 + o(i)) + Etjp?"'' + 0(1)) 

which yields that 

Pi{v) 1, pi{v) ^ 0, i = 2, . . . ,n, as u — > 0. 
In the case oi p = Oj = (c2 — ci)/(a2 — ai) it follows that 

j);c*(l + 0(1)) +p»C''"(l + 0(1)) + E?=3P?!>''C*(1 + 0(1)) 

which gives the limit 

and a similar expression for P2{v). Recall that in this case the coefficient C is found as the 
non-zero solution of 

$(z) = - a,)z^i +p^(62 - aj)z''^ = 0, 

which finally implies that 

, , a,- — 62 / N ttj — bi , , 
^'iW = — P-^ii) = -r — JT^ Pi{t) = 0, ast^oo, 
Oi — 02 02 — Ol 

independently of the initial conditions p^. 

We have that the w-limit set of the replicator equation with the matrix ()4.ip , satisfying the 
genericity condition (|4.5p . consists of the globally attracting equilibrium. This equilibrium is 
either a vertex of the simplex, and in this case the orbit of ()4.8p passing through (1,1) has a 
power asymptote with the exponent given by an index of the second type of the corresponding 
Newton diagram, or this equilibrium is on the 1-skeleton of the simplex, and in this case the 
orbit of (j4.8p passing through (1, 1) has a power asymptote with the exponent given by an index 
of the first type of the corresponding Newton diagram: 

Theorem 2. System (jl.6p with the matrix A given by (j4.ip . satisfying genericity condition 
()4.5p . and such that the conditions of Theorem {1\ hold, always has a global attractor, and this 
attractor can be only an equilibrium. 
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Remark 1. Acc ording to Theorem [T] we excluded the cases when f3j = aj and /3j+i = aj. It 
can be shown (see iBerezovskaval (jl976l )) that in this case we have that system (j4.8|) can have 
asymptotes of the form 

V = Cv'^nnvil + o{l)), 



if — , 



an, and 



Inv 



1 + 0(1)), 



if Pj+i = ctj. These changes do not influence the limits of ()4.9p and hence the conclusion of 
Theorem [2] still holds. 



Remark 2. In the original study lAdams and Sornborger particular ordering of the 

elements of a and b was used: < oi < 02 < • • • < On and 61 > 62 > • • • &n > 0. For the proof 
given above we do not require any particular ordering of the elements of the vectors a and b. 



5 Conclusions 

The replicator equation appears in different problems of evolutionary dynamics in biology and 
economics; it describes the temporary dynamics of frequencies (or probabilities) in heterogeneous 
systems under selective force of natural selection, when the fitness itself is frequency-depended. 
Typically, only limit sets of the systems are under consideration, in particular the rest points 
and their characteristics are an usual object to study. The temporary dynamics of frequencies is 
also of interest and in some applications may be of primary importance. However the problem of 
studying time-dependent behavior is significantly harder than analysis of the limit sets, especially 
for systems of high dimension. 

In this paper we have presented novel methods to analyze the replicator equation (|1.6p . These 
methods, which potentially can be applied to systems of an arbitrary dimension, are based on the 
analysis of the correspon ding selection system, which should be in particular form (|2.ip (for more 
details see Karev ( 20ld )). In Section [2] we provide an algorithmic approach to find the solution 



to the selection system assuming that the initial conditions are given. Our approach consists in 
writing down the corresponding escort system of ordinary differential equations, which in some 
particular cases can be of significantly smaller dimension then the original one. For instance, in 
Example [1] n-dimensional system is replaced with one ordinary differential equation. 

It is worth pointing out that the suggested approach, in addition to the explicit temporal 
dynamics, can be used to infer cj-limit set of the original dynamical system; therefore, we 
only look for attracting asymptotic states and cannot find, e.g., all possible equilibria of the 
replicator equation. In any respect, w-limit sets are what usually is of paramount importance in 
applications because only w-limit sets are what can be observed from the applied point of view. 

For the replicator equation to be suitable for the suggested methods it is usually necessary to 
apply matrix decompositions briefly described in Section [3l One of the possible approaches is the 
singular value decomposition (SVD), which was successfully applied in various static problems to 
reduce the dimension of the data; here we suggest to use SVD for dynamical problems. Generally, 
the escort system, whose asymptotic behavior is of particular interest, is /c-dimensional if the 
original problem has the matrix of rank k. It is therefore tempting to consider an approach 
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when only several singular values are retained, so that we approximate the original matrix of 
rank k with a matrix of rank q < k, for which the escort system if g-dimensional. 

To illustrate the suggested technique we consider a general replicator equation with the 
matrix of rank 1 (Section l3.3|) . Two major finding are that the a;-limit set is always an equilibrium 
and that the existence of the polymorphic (non-isolated) equilibria is not an exception for a 
general matrix A of rank 1. 

As an example of the replicator equation with t he interaction matrix of rank 2 we con- 



sider the problem from lAdams and Sornborgei (120071') . In ge r ieral w e show, using the proposed 

Berezovskava (Il97fil l: iBerezovskava and Kreitseri 



metho ds and the methods of Newton diagram ^ ^ 

(Il976l ). that, for arbitrary dimension and under some suitable conditions (see Theorems [1] and 
[2]), generically one globally stable equilibrium exists on the 1-skeleton of the simplex. We note 
that our conclusions are based on the studying the limit behavior of the solutions of the repli- 
cator equation, which are given in the explicit form, therefore together with the asymptotic 
behavior the time dependent behavior can be effectively analyzed. It is the next step to apply 
the presented methods to the replicator equations with the general interaction matrices of rank 
2 and to higher dimensional problems. 



A Appendix 

Here we show that, although the main theorems for the asymptote s of the trajectories of the 



vector field on a plane were proved only for polynomial systems (see IBerezovskava and Kreitser 



it is straightforward to extend all the results to the system with the right hand sides 



given by quasipolynomials with rational powers. The result follows from the fact that under 
some changes of the variables the indexes of the first and second types of the Newton diagram 
do not change. 

Consider the system 



u = u y piU 'V % 

(A.l) 



where pi, qi G M, Ci,di € Q+. We are interested in power asymptotes u = Cv^il + o(l)) of the 
isolated singular point of this system. Let us make the change of the variables: 

u = y™, V = , m, r € N. 

System (jA.ip takes the form 



r 



^ — 



(A.2) 

X 



It is always possible to choose m and r such that the numbers mdi, rci belong to N, hence 
we can apply the technic of Newton's diagram to system (|A.2p to find the exponents of the 



20 



power asymptotes y = C x ^jl + o(l)). The exponents can be only (see the main text and 
Berezovskava and Kreitsei ( 19761 )) of the form 



r (Ci^-+i - Q.) a _ P'h 



. _ _ V + 1 ^ a. 

Returning to the original variables we obtain that the exponents of the power asymptotes u 
CvP{l + 0(1)) are the indexes of the first or second type 



of the corresponding Newton's diagram built with the rational coordinates dj, Cj. 
Hence the claim is proved. 

Note, that in general we do not need to have identical powers in both equations. 
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